clear;
clc;
tstart=1; %start time
tmax=1000; %max time of simulation
time=tstart:tmax; %time vector

%nparticles = 1000;
nparticles = 1;
npos=2*length(time)+3; %size of pos arrays = 2t+3
pos=1:npos;
nmid_rel=tmax+2; %relative mid position in the below a arrays which corresponds to zero position in the absolute sense
nstart_rel=0;%relative start position in the below a arrays which corresponds to least -ve position in the absolute sense
nend_rel=0;%relative end position in the below a arrays which corresponds to max +ve position in the absolute sense

a_n_u_prev=zeros(1,npos); %a's for spin up at each N at previous time
a_n_d_prev=zeros(1,npos);%a's for spin down at each N at previous time
a_n_u_cur=zeros(1,npos);%a's for spin up at each N at current time
a_n_d_cur=zeros(1,npos);%a's for spin down at each N at current time

%sigma=zeros(1,npos); %the standard deviation of the spatial distribution
sigma=time;
sigma_count = 1;
% T=tmax/2:100:tmax;
T=time;
pos_counter=0;
a=0;

%p_n_tmax_by_2=zeros(1,npos);
alpha=0.025;
% p_n_tmax = zeros(1,npos);
p_n_tmax=zeros(1,npos); %Stores the probability of all positions at a given time t

nGaussian_rand = 10;
% sigma_arr=zeros(1,nGaussian_rand)+sigma_val; %standard deviation of the gaussian random numbers.
% avg=zeros(1,nGaussian_rand); %average for the gaussian random numbers
p_n=zeros(1,nGaussian_rand); %store all the p_n_tmaxs generated for all the gaussian random numbers.

%for all particles
for p=1:nparticles
    %form the initial a vectors for time = 0. consider all the particles to be
    %with spin up
    a_n_u_prev(nmid_rel)=a;
    a_n_d_prev(nmid_rel)=1-a;
    
    %for all times
    for t=time
        display(t);
        pos_counter = 0;
        nstart_rel = nmid_rel - t;
        nend_rel = nmid_rel + t;
        
            %for all possible positions for this particular time, fill up the a
            %arrays for the current time

            for pos_counter=nstart_rel:nend_rel
                        %perform the below calc for a range of random gaussian values
                for p=1:nGaussian_rand
                    a=normrnd(0,alpha);

                    a_n_u_cur(pos_counter)=  ( sqrt(1-a*a) * a_n_d_prev(pos_counter-1) )+ ( a * a_n_u_prev(pos_counter-1) );
                    a_n_d_cur(pos_counter)=( sqrt(1-a*a ) * a_n_u_prev(pos_counter+1)) - ( a * a_n_d_prev(pos_counter+1) );

                      %Below is for the normal quantum walk of the paper without noise.  
%                     a_n_u_cur(pos_counter) =  (1/sqrt(2))*( a_n_d_prev(pos_counter-1) +  a_n_u_prev(pos_counter-1) );
%                     a_n_d_cur(pos_counter)= (1/sqrt(2))*( a_n_u_prev(pos_counter+1) - a_n_d_prev(pos_counter+1) );

                    p_n(p) = a_n_u_cur(pos_counter)^2 + a_n_d_cur(pos_counter)^2; 
        
                end;
                p_n_tmax(pos_counter) = mean(p_n);
        end
        
        %calculate the sigma 
%         if(t==tmax)
%         if(t >= tmax/2 && mod(t,100)==0)
%             sigma = sqrt( (pos^2 .* p_n_tmax) - (pos.*p_n_tmax)^2 )
            sigma(sigma_count) = sqrt( sum(pos.^2 .* p_n_tmax) - (sum(pos.*p_n_tmax)).^2 );
            sigma_count = sigma_count + 1;
%         end;
        
        %set the a arrays of prev time with values of that of the current
        %time and a arrays of current time with zero values
        if(t <tmax)
            a_n_u_prev = a_n_u_cur;
            a_n_d_prev = a_n_d_cur;
            a_n_u_cur = zeros(1,npos);
            a_n_d_cur = zeros(1,npos);
        end;
        
    end;
   
end;

figure(1)
plot(pos-nmid_rel,p_n_tmax);
grid on;

figure(2)
plot(T,sigma);
grid on;

